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Abstract: 

Insulin  is  a  hormone  which  attaches  to  cell  receptors,  allowing  glucose  to  leave  the  blood  stream 
and  provide  energy  to  cells.  Type  I  diabetes  results  when  the  pancreas  does  not  produce  enough 
insulin  to  assist  the  uptake  of  glucose  into  cells.  In  the  case  of  type  II  diabetes,  cell  receptors 
have  decreased  insulin  sensitivity  and  cannot  use  insulin  to  transport  glucose  into  cells  properly. 
Both  conditions  present  the  danger  of  causing  unhealthy  glucose  levels  in  the  blood  stream  and 
are  treated  with  insulin  injection  therapy  to  trigger  glucose  uptake  in  the  cells. 

A  mathematical  model  of  natural  glucose  and  insulin  control  in  the  body  allows  for  a  deeper 
understanding  of  the  malfunctions  that  cause  diabetes,  as  well  as  a  more  quantified  approach  to 
regulating  blood  glucose  levels  with  insulin  injection  therapy.  Cobelli  et  al.  presented  a 
simulation  model  which  describes  the  kinetics  of  glucose  digestion  and  absorption  that  occur 
after  a  meal.1 

First,  this  research  seeks  to  extend  the  Cobelli  model  of  glucose  and  insulin  dynamics  to  include 
both  long  and  short-acting  insulin  inputs  currently  used  to  treat  diabetic  patients.  Secondly,  the 
research  will  introduce  a  personalized  approach  to  treatment  by  adapting  the  combined  model 
parameters  over  time  in  response  to  observed  patient  feedback  data.  As  a  result,  the  updated 
Cobelli  model  parameters  can  be  fitted  to  the  individual  patient.  The  resulting  model  can  be 
used  in  future  research  to  predict  the  outcomes  of  different  insulin  dosages  and  determine  the 
best  treatment  option  to  achieve  desired  glucose  levels  in  diabetic  patients.  Consequently  this 
research  forms  the  foundation  for  a  personalized  prescription  tool  to  improve  diabetic  patient 
care. 

Keywords:  insulin,  glucose,  diabetes,  injection,  model 


Enclosures:  (1)  Cobelli  Model  Flow  Chart 

(2)  Combined  Cobelli  Model  Equation  Layout 

(3)  Cobelli  Parameter  Spreadsheet 

(4)  Randomized  Patient  Data  Matlab  Code 

(5)  Parameter  Sensitivity  Testing  Matlab  Code 

(6)  Parameter  Adaptation  Matlab  Code 
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Understanding  Insulin  and  Diabetes 

Human  cells  use  insulin  to  extract  glucose  from  the  bloodstream,  using  the  simple  sugar  for 
energy.  The  hormone  insulin  is  analogous  to  a  key  that  is  needed  to  “unlock”  cells  by  attaching 
to  cell  receptors,  allowing  glucose  to  enter.  Malfunctions  in  the  regulation  of  insulin  and  glucose 
in  the  body  can  give  rise  to  pathological  conditions  including  both  type  I  and  type  II  diabetes. 
Type  I  diabetes  results  when  the  pancreas  does  not  produce  enough  insulin  to  assist  the  uptake  of 
glucose  into  cells.  Type  II  diabetes  is  a  condition  in  which  cell  receptors  have  decreased  insulin 
sensitivity  and  cannot  use  insulin  to  intake  glucose  properly.  Both  conditions  result  in  abnormal 
glucose  levels  in  the  blood  stream,  potentially  resulting  in  complications  including  heart  disease, 
stroke,  high  blood  pressure,  blindness,  kidney  failure,  and  nervous  symptom  damage. 

Both  type  I  and  type  II  diabetes  are  currently  treated  with  insulin  injection  therapy  to  trigger 
glucose  uptake  in  the  cells.  The  Juvenile  Diabetes  Research  Foundation  launched  the  Artificial 
Pancreas  Project  in  2006  to  push  for  the  development  of  an  artificial  pancreas. 1V  This  implant 
would  work  in  the  place  of  a  type  1  diabetic’s  failing  pancreas,  regulating  blood  glucose  levels  in 
the  blood  by  administering  the  appropriate  dose  of  insulin. 

Project  Goals 

This  project  will  bridge  the  gap  between  current  artificial  pancreas  research  and  clinical 
application  to  improve  treatment  methods  for  both  type  I  and  type  II  diabetic  patients.  While  the 
development  of  an  artificial  pancreas  is  a  desirable  long-term  goal,  a  physiologically  accurate 
glucose  control  model  has  immediate  potential  as  a  prescription  tool  for  physicians  to  use.  This 
prescription  tool  will  help  physicians  to  better  understand  their  specific  patient’s  condition  and 
determine  the  best  course  of  insulin  therapy  needed  to  achieve  healthy  blood  glucose  levels. 

Figure  1  outlines  the  research  objectives  in  a  flow  chart  format.  First,  the  research  will  improve 
the  clinical  relevance  of  an  existing  model  of  glucose  dynamics  which  currently  only  accounts 
for  the  influence  of  insulin  that  is  naturally  secreted  from  the  pancreas.  The  improved  model  will 
include  additional  insulin  inputs  currently  used  in  injection  therapy.  The  second  objective  is  to 
adapt  the  Cobelli  model  parameters  to  investigate  a  personalized  approach  to  treatment,  fine 


3 


tuning  the  model  to  best  represent  the  patient’s  condition.  These  two  research  steps  will  prepare 
the  combined  Cobelli  and  Insulin  models  for  future  use  in  adaptive  model  predictive  control, 
building  the  foundation  for  a  personalized  prescription  tool  for  physicians  to  use.  The  Juvenile 
Diabetes  Research  Foundation  states  the  “development  of  novel  leaming/patient  specific 
artificial  pancreas  control  algorithms  in  closed  loop  control  systems”  as  one  of  its  top  priority 
areas  for  the  fiscal  year  of  2013.vThis  research  will  contribute  to  an  individualized  prescription 
tool  that  physicians  can  use  to  better  treat  diabetic  patients,  while  furthering  the  development  of  a 
control  algorithm  that  could  be  applied  to  a  functioning  artificial  pancreas  in  the  future. 


Objectives  to  Achieve 


Objective  I: 

Expand  model 
to  include 
slow  and  fast 
acting  insulin 
inputs 


Objective  III: 

Objective  II:  \ 

Investigate 

Adapt  the 

insulin  input  J  <— 

model 

prediction 

parameters 

Figure  1:  Research  Objectives 

Research  Timeline 

My  first  semester  research  goals  were  centered  on  developing  a  working  understanding  of  the 
Cobelli  model,  solving  the  Cobelli  differential  equations  using  Simulink,  finding  an  appropriate 
mathematical  model  of  insulin  injections  to  incorporate  into  the  Cobelli  model,  and  investigating 
the  Cobelli  model  parameter  sensitivity  to  prepare  for  parameter  adaption. 

During  second  semester,  a  virtual  patient  was  created  using  the  Cobellii  model  with  a 
randomized  set  of  model  parameters.  This  generated  patient  data  was  used  for  analysis 
throughout  the  research,  as  patient  data  from  Anne  Arundel  Medical  Center  was  obtained  in 
early  April  and  could  not  be  used  in  the  project.  Parameter  sensitivity  analysis  results  identified  a 
subset  of  the  Cobelli  model  parameters  that  would  be  most  effective  in  adapting  the  published 
Cobelli  diabetic  patient  parameters  to  fit  the  virtual  patient  data.  Additionally,  models  of  slow 
and  fast  acting  insulin  injections  were  added  to  the  adaptive  Cobelli  model.  This  adaptive  model 
will  be  useful  in  future  research  to  create  a  personalized  prescription  tool  to  improve  the 
treatment  of  diabetic  patients.  The  remainder  of  this  report  will  describe  the  pertinent 
background  information  and  results  of  my  work  throughout  the  year. 

Cobelli  Background  and  Modeling 

Cobelli  et  al.  presented  a  nonlinear  simulation  model  composed  of  glucose  and  insulin 
subsystems  and  four  unit  process  models  which  describe  the  kinetics  and  physiological  events  of 
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glucose  digestion  and  absorption  in  key  organs  after  a  meal.  They  used  a  "triple  tracer  meal 
protocol"  that  provided  estimates  of  major  glucose  and  insulin  fluxes  (i.e.  flow  rates)  during  the 
course  of  a  meal.vl  The  model  relates  glucose  and  insulin  fluxes  to  compartmental  concentrations 
in  order  to  characterize  the  flow  through  each  subsystem,  using  28  differential  equations  to 
represent  glucose  and  insulin  dynamics  in  a  physiologically  accurate  manner.  Cobelli’s  model  is 
an  appealing  foundation  for  my  research  because  a  computer  simulator  of  type  1  diabetes  was 
based  on  this  model  and  approved  by  the  FDA  as  a  substitute  to  animal  trials  for  preclinical  drug 
testing.  vn 

For  the  sake  of  brevity,  only  the  differential  equations  characterizing  the  insulin  subsystem  will 
be  described  further.  The  two-compartment  insulin  subsystem,  displayed  in  Figure  2,  describes 
the  response  of  the  insulin  masses  in  the  liver  (7;  -  Equation  1 )  and  plasma  (Ip  —  Equation  2)  to 
insulin  secreted  by  beta  cells  in  the  pancreas  (S).  Specifically,  these  masses  vary  as  insulin  is 
utilized  in  the  periphery  (muscle  and  tissue)  cells  to  assist  in  glucose  uptake  or  is  lost  within  the 
liver.  The  two  equations  shown  below  describe  the  flow  (change  in  mass  over  time)  of  insulin 
through  each  compartment.  Note  that  the  dot  notation  indicates  the  derivative  with  respect  to 
time.  The  variables  ,  m2 ,  m, ,  and  w4  are  rate  parameters  describing  the  exchange  between  the 
liver  and  plasma  insulin  compartments  as  well  as  the  rate  at  which  insulin  is  either  lost  within  the 
liver  or  used  in  surrounding  cells. 


Figure  2:  Insulin  Subsystem  and  Equations  from  Cobelli  Model 

h( 0  =  S(t)  +  (m2Ip(t)  -  mj/^t))  —  m3/;(t)  ( liver  insulin )  Equation  1 

Ip(t)  =  (m^iit)  —  m2/p(t))  —  tn4Ip(t)  ( plasma  insulin )  Equation  2 

Before  proceeding,  I  wanted  to  possess  an  in-depth  understanding  of  Cobelli’s  glucose  and 
insulin  model.  To  do  so,  I  created  three  reference  documents  including  a  block  diagram  relating 
the  Cobelli  model  subsystems  (Enclosure  1),  a  similarly  structured  diagram  outlining  the 
necessary  equations  from  Cobelli’s  original  paper  (Enclosure  2),  and  a  chart  of  all  model 
parameters  and  their  physiological  meaning  (Enclosure  3).  MIDN  1/C  Brandon  Meek,  a  systems 
engineering  honors  student  working  on  a  related  capstone  project,  then  used  these  reference 
documents  to  build  a  Simulink  model  of  the  Cobelli  equations.  Simulink  is  a  computer  program 
available  at  USNA  that  uses  a  block  diagram  (graphical)  environment  to  solve  differential 
equations  and  simulate  dynamic  systems  over  time.  Figure  3  shows  the  Simulink  model 
structure.  Each  block  in  this  image  is  a  subsystem  containing  the  relevant  subsystem  equations. 
The  insulin  subsystem  is  circled. 
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Figure  3:  Cobelli  Simulink  Model 


The  interior  of  the  insulin  subsystem  is  shown  in  Figure  4.  The  two  circled  integrator, 
blocks  indicate  that  two  differential  equations  are  solved  within  the  subsystem. 


Figure  4:  Simulink  Insulin  Subsystem 

The  scope  allows  the  user  to  plot  the  plasma  insulin  mass  over  time  and  compare  the  results  to 
those  published  in  Cobelli ’s  paper/111  The  graph  in  Figure  5  represents  the  Simulink  insulin  data 
as  seen  through  the  MATLAB  scope  after  a  patient  is  given  a  meal.  As  expected,  the  insulin 
level  increases  to  help  transport  the  ingested  glucose  from  the  blood  stream  into  cells  to  be  used 
for  energy.  This  graph  matches  the  insulin  data  published  by  Cobelli,  indicating  that  the  Cobelli 
model  can  be  accurately  simulated  and  manipulated  using  Simulink  code. 
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Figure  5:  Simulink  Insulin  Data 

Once  the  Simulink  Cobelli  model  was  verified  to  match  the  published  results,  it  was  used  to 
create  a  “virtual  patient”  by  randomizing  the  model  parameters.  Cobelli  et  al.  published  model 
parameters  for  both  healthy  and  diabetic  patients. lx  The  percent  change  between  these  values 
was  calculated,  and  the  parameters  were  randomized  within  a  range  of  plus  or  minus  10%  of  the 
diabetic  patient  parameter  value.  Using  this  interval,  multiple  unique  patient  parameter  sets  could 
be  used  to  generate  data  using  the  Cobelli  Simulink  model.  Figure  7  shows  the  simulation  results 
for  the  one  meal  using  the  original  Cobelli  parameters,  as  well  as  three  different  patients  with 
randomized  parameter  values. 
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Randomized  Virtual  Patient  Glucose  Data 


Figure  5:  Virtual  Patient  Data 

Objective  1:  Improve  Clinical  Relevancy  by  Expanding  the  Cobelli  Model 

Diabetes  patients  are  currently  treated  with  both  rapid- acting  and  long-acting  insulin  injections  to 
achieve  effective  glucose  control  throughout  the  course  of  a  day.  The  characteristics  of  common 
insulin  analog  injections  are  displayed  below.x 


Rapid  Acting 

Onset 

Peak 

Duration 

Humalog  or  Lispro 

15-30  min 

30-90  min 

3-5  hr 

Novolog  or  Aspart 

10-20  min 

45-50  min 

3-5  hr 

Apidra  or 

Glulisine 

20-30  min 

30-90  min 

1-2.5  hr 

Long  Acting 

Onset 

Peak 

Duration 

Lantus 

1-1.5  hr 

None  (steady) 

20-24  hr 

Levemir 

1-2  hr 

6-8  hr 

Up  to  24  hr 

Table  1:  Common  Insulin  Injection  Profi 

les 
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The  Cobelli  model  must  be  modified  to  include  both  rapid-acting  and  long-acting  insulin  inputs 
to  the  plasma  insulin  compartment.  This  modification  will  make  the  Cobelli  model  more 
relevant  to  doctors  considering  the  model  as  a  potential  prescriptive  aid.  The  altered  insulin 
subsystem  is  shown  in  Figure  7. 


Figure  6:  Insulin  Subsystem  with  added  injections 


The  original  research  proposal  suggested  the  use  of  first  order  transfer  functionswith  differing 
time  constants  and  an  appropriate  time  delay  for  rapid-  and  long-acting  insulin  injections,  as 
shown  in  Equation  3  and  4. 


Slow  Acting  Insulin  Transfer  Function 


400 

In(s) 


1 

TsS  + 


Fast  Acting  Insulin  Transfer  Function 


his o 

7n(s) 


T/ 


—  ■  (e  tdS )  Equation  3 

1 

—  -  Equation  4 


However,  further  literature  study  revealed  significant  prior  work  on  the  subject  of  insulin  analog 
injection  modeling.  The  literature  trail  showed  that  a  single  first  order  equation  oversimplifies 
the  injection  and  absorption  of  insulin  under  the  skin,  through  the  capillary  membrane,  and  into 
the  bloodstream.  A  review  of  four  major  papers  resulted  in  the  final  representation  chosen.  My 
findings  are  summarized  below. 


Similar  to  our  proposal,  Huang  et  al.  presents  a  simple  first  order  model  to  describe  the  dynamics 
of  impusive  insulin  injections.5"  However,  this  model  does  not  have  separate  compartments  for 
the  tissue  and  blood  stream  with  associated  time  delays.  Next,  a  three  compartment  model 
presented  by  Shimoda  et.  al.  was  considered.5"1  This  model  includes  a  separate  injection  depot 
compartment,  which  is  represented  with  a  first  order  delay.  However,  the  model  only  applies  to 
rapid  acting  insulin,  as  a  first  order  delay  is  not  enough  to  simulate  the  slow  dissolving  process  of 
long  acting  injections. 


The  most  helpful  paper  throughout  this  literary  search  was  a  review  of  multiple  insulin  injection 
models  written  by  G.  Nucci  and  C.  Cobelli.51111  In  this  paper,  Cobelli  dismisses  the  Shimoda 
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model  because  it  was  conceived  in  the  1980’s,  before  modem  insulin  analogs  were  available.  A 
more  promising  model  presented  by  Trajanoski  et.  al.  describes  the  breakdown  of  insulin  analogs 
from  an  original  hexameric  state  to  a  dimeric  state  that  can  be  absorbed  into  the  bloodstream. 

The  model  includes  an  imaginary  “bound  state”  to  represent  the  delayed  dynamics  of  long  acting 
insulin  injections.  Cobelli  praises  this  model  for  its  chemical  accuracy,  dependency  on  both  time 
and  injection  depth,  and  inclusion  of  an  inverse  relationship  between  injection  volume  and 
absorption  rate.  At  the  end  of  the  review  paper,  Cobelli  concludes  by  suggesting  a  simplification 
of  the  Trajanoski  model  structure  to  reduce  computational  burden. 

Jiaxu  Li  and  Yang  Kuang  answered  with  their  simplified  version  of  Trajanoski’ s  model.xlv  The 
model  equations  vary  with  time  (rather  than  both  time  and  injection  depth),  include  a 
biologically  accurate  Michaelis-Menten  description  of  insulin  degradation,  and  lessen  the 
computational  burden  of  the  original  model.  This  model  appears  to  be  the  best  fit  to  incorporate 
into  Cobelli’ s  model  of  insulin  and  glucose  dynamics:  it  describes  the  chemistry  of  insulin 
injection  absorption,  appropriately  accounts  for  compartmental  delays,  and  has  been  reviewed 
positively  by  Cobelli.  The  rapid  (short)  acting  and  long  acting  injection  model  equations  shown 
below  were  modeled  in  Simulink  to  be  incorporated  into  the  Cobelli  model  insulin  equation.  In 
these  equations,  Itotai  represents  the  total  plasma  insulin  mass,  while  fast  and  Isiow  represent  the 
individual  insulin  injection  concentrations.  xv 


Short-acting  insulin:  Lispro 

Hfast (0  =  —  p{H(t)  —  qD3(t ))  ( hexameric  form )  Equation  5 


✓  ,  x  bD(t ) 

DfastW  =  p(H(t)  -  qD3(t))-——K 

J-Ti/; 


(i dimeric  form )  Equation  6 


rbDit ) 


' total  (0 

lfast(t)  =  ~  dilfast(t )  ( plasma  concentration  of  fast  insulin)  Equation  7 

Hfasti 0)  =  injection  amount;  Dfast( 0)  =  0;Ifast(0)  =  lfast 0 


Long-acting  insulin:  Glargine 


B'(t)  =  —kB(t)  * 


1  +  Hit) 


(bound  form)  Equation  8 


Hsiow(.t)  =  —p(.H(t)  ~  qD3(t ))  +  kB(t)  *- — (hexameric  form)  Equation 9 

1  +  H(t) 


bD(t) 


Dsiow(t)  =  p(//(t)  —  qD3(t ))  —  - - - —  (dimeric  form)  Equaution  10 

1  +  hotai(t ) 


(t) 


flowif) 


rbD{t) 

1+hotaiO 


5(0)  = 


dilfastit )  (plasma  concentration  of  fast  insulin) 
injection  amount ;  Dsiow( 0)  =  0;  Isiow( 0)  =  Isiow 0 


Equation  11 
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Equations  7  and  1 1  represent  the  plasma  insulin  concentration  of  fast  and  slow  acting  insulin 
which  Li  and  Kuang  proposed.  When  the  fast  and  slow  acting  insulin  models  were  originally 
added  to  the  Cobelli  plasma  insulin  equation  (represented  in  Equation  2),  the  resulting  equation 
was  as  follows: 

l total (0  =  (mi It(t)  -  m2/p(t))  -  m4/p(t)  (plasma  insulin )  +  Ifast(t) 

+  Isiow(t )  Equation  12 

However,  there  are  two  key  differences  between  the  Cobelli  and  Li/Kuang  models  that  are 
important  to  note.  First,  the  insulin  equations  are  represented  using  different  units,  uU/ml  in 
Li/Kuang  and  pmol/kg  in  Cobelli.  Additionally,  this  difference  in  units  reflects  the  fact  that  the 
insulin  injection  equations  represent  insulin  concentrations  in  the  bloodstream,  while  the  Cobelli 
equations  describe  a  mass  of  plasma  insulin.  To  remedy  these  differences  and  accurately 
combine  the  Cobelli  and  insulin  injection  models,  the  combined  insulin  equation  was  altered  to 
multiply  each  insulin  injection  term  by  the  total  plasma  insulin  volume  in  the  Cobelli  model  (Vi). 
This  correction,  along  with  additional  unit  conversions  within  the  Li/Kuang  insulin  model 
parameters,  placed  the  combined  equation  terms  consistently  in  pmol/kg  to  represent  the 
combined  plasma  insulin  mass. 

A  few  additional  changes  were  made  to  the  Li/Kuang  insulin  injection  equations  to  form  an 
appropriate  total  plasma  insulin  equation.  Originally,  the  proposed  equations  for  the  fast  and 
slow  acting  insulin  injections  (Equations  5  and  9)  were  dependent  on  the  total  plasma  insulin 
concentration,  both  in  the  insulin  addition  term  and  the  negative  insulin  degradation  term.  The 
equations  were  edited  so  that  the  addition  of  dimeric  insulin  is  inversely  proportional  to  the  total 
insulin  mass,  ltotai(X)^  but  the  degradation  of  insulin  is  dependent  on  the  concentration  of  the 
injection,  //ast(0-  The  initial  condition  of  the  slow  and  fast  acting  insulin  concentration 
equations  (Equations  7  and  1 1  -//ast0  and  /s;ow0)  was  originally  set  at  6  uU/ml  to  normalize  the 
patient’s  plasma  glucose  level  to  basal  conditions.  This  was  because  patients  were  given  insulin 
overnight  to  normalize  their  plasma  glucose  levels,  isolating  the  effect  of  the  insulin  injections 
and  providing  data  that  would  prove  useful  in  modeling  the  injection  dynamics. xvlxvn  When  the 
insulin  injection  concentrations  were  added  to  the  total  insulin  mass  in  Equation  10,  the  insulin 
initial  conditions  were  set  to  zero.  The  insulin  component  of  the  Cobelli  model  will  account  for 
any  plasma  insulin  produced  by  the  body,  so  no  baseline  initial  condition  is  needed. 

Itotaii 0  =  (mi/jC t)  -  m2/p(t))  -  m4/p(t)  ( plasma  insulin )  +  7/ast(t)  *  Vt  +  /siow(t)  *  vi 

Equation  10 

The  combined  model  was  coded  in  Simulink  and  used  to  administer  both  slow  and  fast  acting 
insulin  to  a  virtual  patient  over  a  24-hour  period.  The  patient  received  40  units  of  slow  acting 
insulin  at  time  zero,  along  with  his  first  meal.  He  then  received  8  units  of  fast  acting  insulin 
along  with  his  second  meal  6  hours  later.  Here,  1  unit  of  insulin  equals  6.94x  106  pmol/1.  Figure 
9  shows  the  resulting  glucose  levels  of  the  patient,  as  compared  to  the  same  patient  who  received 
no  insulin  or  only  the  slow  insulin  injection  at  time  zero.  When  slow  insulin  is  administered,  it 
lowers  the  patient’s  peak  and  resting  glucose  levels  as  shown  by  the  red  line  graph.  The  fast 
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insulin,  administered  along  with  meal  two,  creates  an  additional  decrease  in  the  blood  glucose 
level  that  could  not  be  achieved  with  slow  insulin  alone.  This  is  best  seen  by  comparing  the  red 
(slow  and  fast  insulin)  and  green  (slow  insulin  only)  lines  in  in  the  second  spike  of  glucose 
concentration.  In  this  case,  the  patient  was  given  five  times  the  amount  of  slow  acting  insulin  as 
compared  to  fast  insulin,  so  it  is  not  surprising  that  the  fast  acting  insulin  does  not  have  a  large 
impact  on  the  patient’s  glucose  levels. 

Glucose  Concentration:  Slow  insulin  =  40  units,  Fast  insulin  =  8  units 


Time  (min) 


Figure  9:  Combining  Cobelli  and  Insulin  Models 

Figure  10  shows  the  insulin  concentration  for  the  same  patient,  as  slow  insulin  is  administered  at 
time  0  and  fast  insulin  is  administered  six  hours  later.  The  injected  insulin  is  necessary  to  help 
the  ingested  glucose  leave  the  blood  stream  and  enter  the  cells  following  a  meal.  Focusing  on  the 
valleys  following  each  meal,  the  insulin  injections  prevent  the  patient’s  insulin  levels  from 
dropping  as  low  as  they  would  without  injected  insulin.  These  increased  insulin  levels  assist  in 
glucose  uptake,  lowering  plasma  glucose  levels  between  meals. 
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Insulin  Concentration:  Slow  insulin  =  40  units,  Fast  insulin  =  8  units 


Figure  10:  Showing  Insulin  Concentration  changes 

This  combined  model  seems  to  accurately  represent  the  desired  effect  of  injecting  a  patient  with 
external  insulin:  lowering  blood  glucose  levels,  particularly  between  meals.  These  results  will 
need  to  be  compared  to  patient  glucose  readings  over  the  course  of  insulin  treatment  to  verify  the 
accuracy  of  the  combined  model. 


Objective  2:  Adapt  the  Cobelli  Model  Parameters 

Once  the  Cobelli  model  had  been  verified  with  Simulink  and  combined  with  an  appropriate 
insulin  injection  model,  my  focus  shifted  to  preparing  to  adapt  the  original  Cobelli  model 
parameters  to  best  fit  the  condition  of  an  individual  patient.  The  first  step  was  to  evaluate  the 
sensitivity  of  the  glucose  and  insulin  responses  of  the  Cobelli  model  to  parameter  changes.  This 
analysis  provided  a  subset  of  the  most  sensitive  model  parameters  that  would  be  most  effective 
in  parameter  adaptation  to  create  a  personalized  Cobelli  model.  The  MATLAB  code  in  Enclosure 
4  changes  each  of  the  35  model  parameters  by  incremental  percentages  of  their  original  value, 
runs  the  Cobelli  Simulink  model,  and  obtains  the  resulting  glucose  data.  Then,  the  code 
calculates  the  error  between  the  altered  and  original  response  to  compute  the  model’s  sensitivity 
to  the  parameter  change.  The  average  sensitivity  was  determined  for  each  parameter  using 
equation  1 1 . 


Sens  =  100  * 


N 


( altered  —  original )2 
#  readings 


ma x(gluc)  —  min  (glue)  parameter  change 


Equation  11 


13 


Several  different  methods  of  parameter  sensitivity  testing  were  used  to  ensure  the  proper 
parameters  were  chosen  for  adaptation.  The  sensitivity  was  first  calculated  using  error  in 
resulting  glucose  data,  followed  by  the  insulin  data  error,  as  well  as  a  sum  squared  error 
combining  both  glucose  and  insulin  results.  Four  parameters  (Fens,  gamma,  Heb,  and  f)  did  not 
change  between  healthy  and  diabetic  patients  in  Cobelli’s  paper,  so  these  parameters  will  not  be 
adapted  to  fit  patient  data  and  have  been  removed  from  all  data  tables  below.  Table  2  shows  the 
parameters  which  were  most  sensitive  when  both  insulin  and  glucose  error  were  considered. 


Parameter 

Name 

Parameter  Description 

Average 

Sensitivity 

kpl 

EGP  at  zero  glucose  and 
insulin 

10711.849 

Vg 

Distribution  volume  of 
glucose 

1538.5491 

kp4 

Amplitude  of  portal 
insulin  action  on  liver 

888.2471 

kl 

Glucose  rate  parameter 
(b/n  plasma  and  tissue) 

677.07393 

Vi 

Distribution  volume  of 
insulin 

474.51345 

m6 

Rate  parameter  for 
Hepatic  extraction 

382.68959 

m4 

Insulin  rate  parameter 
(periphery  degredation) 

284.53814 

Table  2:  Parameter  Sensitivity  Results  -  Combining  Insulin  and  Glucose  Error 

Since  the  parameter  sensitivity  is  determined  using  glucose  and  insulin  data,  it  is  logical  for  the 
most  sensitive  parameters  to  be  involved  with  endogenous  glucose  production  from  the  liver 
(kpl  and  kp4),  parameters  within  the  glucose  subsystem  (Yg,  kl),  or  within  the  insulin 
subsystem  (Yi,  m4,  and  m6). 

After  this  sensitivity  testing  was  conducted,  additional  sensitivity  testing  was  designed  to  not 
only  increase  the  parameter  values  and  determine  resulting  sensitivity,  but  to  decrease  the 
parameter  values  as  well.  This  sensitivity  testing  yielded  the  results  shown  in  Table  3.  These 
results  have  overlapping  sensitive  parameters  with  those  in  Table  2,  except  for  the  three 
parameters  highlighted  in  blue. 


Parameter 

Name 

Parameter  Description 

Average 

Sensitivity 

kpl 

EGP  at  zero  glucose  and 
insulin 

5777.668044 

m5 

Rate  parameter  for 
Hepatic  extraction 

1722.260222 
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kabs 

Rate  constant  of  intestinal 
absorption 

977.5885681 

kgri 

Rate  of  grinding 

795.9737116 

Vg 

Distribution  volume  of 
glucose 

533.1015597 

kp4 

Amplitude  of  portal 
insulin  action  on  liver 

474.7515014 

kl 

Glucose  rate  parameter 
(b/n  plasma  and  tissue) 

416.7587145 

m6 

Rate  parameter  for 
Hepatic  extraction 

344.0057034 

Table  3:  Parameter  Sensitivity  Results  -  Using  +  and  -  parameter  changes 

After  this  first  two  rounds  of  parameter  sensitivity  testing,  it  was  decided  that  parameter 
adaptation  would  be  evaluated  first  by  adapting  the  top  five  parameters  from  the  combined 
insulin  and  glucose  sensitivity  testing  (kpl,  Vg,  kp4,  kl,  VI).  Then,  parameter  adaptation  would 
be  evaluated  when  ten  parameters  are  adapted,  including  the  top  seven  combined  sensitivity 
results  as  well  as  the  three  additional  parameters  highlighted  in  table  3. 

The  parameter  adaptation  code  (Enclosure  5)  is  a  Matlab  script  which  uses  an  optimization 
function  to  choose  the  best  set  of  parameters  to  match  patient  glucose  data.  The  script  runs 
through  multiple  iterations,  trying  different  parameter  sets  within  a  determined  set  of  upper  and 
lower  bounds,  solving  the  Cobelli  model  equations,  and  computing  the  sum  squared  error  (or 
cost)  between  the  simulated  data  and  the  provided  patient  data.  In  this  case,  the  patient  data  was 
generated  using  a  Simulink  model  with  randomized  parameters,  as  described  before  with  the  best 
fit.  Once  the  program  can  no  longer  lower  the  cost,  it  returns  the  best  fit  set  of  parameters.  The 
optimization  function  also  returns  the  cost  between  the  patient  glucose  data  and  the  data 
simulated  using  the  chosen  parameters.  Here,  the  goal  is  to  reduce  the  cost  without  adapting  all 
35  Cobelli  model  parameters,  thus  reducing  computational  burden. 

Figure  1 1  shows  the  results  of  parameter  adaptation  in  which  all  35  Cobelli  model  parameters 
were  adapted  to  fit  the  randomized  patient  data.  In  this  case,  the  glucose  data  simulated  with  the 
adapted  parameters  (black  line)  and  the  patient  data  (blue  line)  completely  overlap.  This  and  the 
extremely  low  calculated  cost  prove  that  adaptation  was  successful.  However,  the  parameter 
adaptation  code  ran  for  over  10  minutes  before  determining  the  optimal  set  of  parameters.  In 
future  research,  parameter  adaptation  will  be  a  portion  of  a  larger  input  prediction  routine,  and  it 
is  unreasonable  to  allow  this  much  time  for  adaptation.  These  results  confirmed  the  need  to 
choose  a  smaller  set  of  parameters  to  more  efficiently  adapt  the  Cobelli  model  to  fit  patient  data. 
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Paramptpr  AHantatinn  RpmiIIw  all  35  Paramptpr<; 


Figure  11:  Adapting  all  35  parameters 

Next,  the  parameter  adaptation  code  was  written  to  adapt  the  top  five  sensitive  parameters  (as 
determined  from  the  combined  insulin  and  glucose  data)  and  hold  all  other  parameters  constant 
at  the  Cobelli  values.  Doing  so  greatly  increased  the  cost  (from  0.217  to  961).  To  put  the 
resulting  cost  in  perspective,  I  compared  the  parameter  adaptation  cost  (961)  to  the  cost  (or  sum 
squared  error)  between  the  patient  data  and  the  Cobelli  model  glucose  response  using  the 
published  diabetic  patient  parameters.  The  adapted  cost  is  still  less  than  8%  of  the  cost  calculated 
using  the  data  generated  from  the  original  Cobelli  parameters;  These  results  are  shown  in  Figure 
12  below. 
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Parameter  Adaptation  Results:  Top  5  Parameters 


Figure  12:  Adapting  top  5  parameters  (glucose  and  insulin  sensitivity) 

The  next  round  of  parameter  fitting  changed  ten  parameters  to  fit  the  randomized  patient  data.  As 
shown  in  Figure  13,  the  cost  (161)  is  significantly  reduced  by  allowing  just  five  more  parameters 
to  adapt  and  fit  the  patient  data.  This  new  cost  is  16.7%  of  the  cost  resulting  from  adapting  only 
the  top  five  parameters.  The  code  ran  in  a  reasonable  amount  of  time,  suggesting  that  this  is  a 
reasonable  number  of  parameters  to  adapt  in  order  to  personalize  the  Cobelli  model. 
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Parameter  Adaptation  Results:  Top  10  Parameters 


Figure  13:  Adapting  top  10  Parameters 

After  these  first  two  rounds  of  parameter  adaption,  we  revisited  the  parameter  sensitivity  testing 
with  a  different  approach  to  see  if  the  adaptation  results  could  be  improved.  This  time,  instead  of 
altering  each  parameter  by  a  percentage  of  its  original  value,  we  based  the  parameter  changes  on 
the  percent  change  between  the  healthy  and  diabetic  value  for  each  parameter.  This  approach 
was  designed  to  emulate  the  method  we  used  to  generate  the  randomized  patient  data,  by  altering 
the  original  Cobelli  diabetic  value  by  a  random  percentage  of  the  percent  change  between 
healthy  and  diabetic  parameter  values.  Each  parameter  was  changed  by  positive  and  negative  5% 
and  10%  of  the  percent  change  between  healthy  and  diabetic  parameter  values.  The  sensitivity 
was  then  calculated  using  the  resulting  glucose  data,  and  the  top  five  and  ten  parameters  were 
chosen  for  additional  rounds  of  adaptation  testing. 

Table  4  summarizes  the  parameter  adaptation  results  for  all  5  scenarios:  adapting  all  35 
parameters,  the  top  5  parameters  from  the  combined  glucose  and  insulin  data,  top  10  parameters 
(including  positive  and  negative  parameter  changes),  the  top  5  “new”  parameters  chosen  by 
considering  the  percent  change  in  parameter  values,  as  well  as  the  top  10  parameters  from  the 
same  round  of  parameter  sensitivity  testing.  The  first  two  columns  compare  the  glucose  cost 
calculated  using  the  adapted  parameter  values  to  the  cost  generated  using  the  original  diabetic 
Cobelli  parameter  values.  The  new  method  of  parameter  sensitivity  testing  picked  out  parameters 
which  caused  an  increased  cost  as  compared  to  the  original  method.  The  original  top  ten 
parameters  are  therefore  the  most  effective  in  adapting  to  personalize  the  Cobelli  model  to  fit 
patient  glucose  data. 
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The  third  and  fourth  columns  of  Table  4  record  the  sum  of  the  percent  error  between  the  adapted 
and  original  patient  parameters,  as  well  as  the  percent  error  between  the  Cobelli  parameters  and 
actual  patient  parameters.  This  is  a  potentially  helpful  analysis  tool  because  it  shows  how  close 
the  parameters  were  to  matching  the  randomized  patient  values.  All  five  cases  have  very  low 
parameter  error  percentages,  making  this  a  difficult  tool  to  use  to  distinguish  between  adaptation 
methods.  Additionally,  the  parameter  percent  error  did  not  always  decrease  when  the  parameters 
were  allowed  to  adapt.  The  most  important  metric  to  use  when  analyzing  parameter  adaptation 
effectiveness  is  the  glucose  cost,  because  patient  glucose  data  will  be  used  to  personalize  the 
Cobelli  model  and  determine  insulin  injections  in  any  prescription  tool.  The  individual  adapted 
parameters  may  not  come  closer  to  the  patient  parameter  values,  but  this  is  of  little  importance  as 
long  as  the  simulated  glucose  data  fits  the  provided  patient  data  with  as  low  of  cost  as  possible. 


Adapted  glucose 

cost 

Diabetic  glucose 

Cost 

%  error  adapted 
parameter 

%error 

diabetic 

parameter 

1  meal,  all  35  parameters 

0.217 

12100 

2.07939675 

1.983223512 

1  meal,  5  parameters 

961 

12100 

1.807557703 

1.828632793 

1  meal,  10  parameters 

119 

12100 

1.426699217 

2.491692 

1  meal,  5  new  parameters 

1630 

12100 

0.686665424 

8.012487486 

1  meal,  10  new  parameters 

1046.7 

12100 

1.520362334 

0.627409478 

Table  4:  Summarizing  Adaptation  Results 
Conclusion  and  Future  Directions 

The  two  main  objectives  of  this  research  project  were  accomplished  over  the  course  of  the  year: 
the  Cobelli  model  was  combined  with  physiologically  accurate  models  of  slow  and  fast  acting 
insulin  injections.  The  parameters  of  the  Cobelli  model  equations  were  adapted  to  personalize  the 
model,  fitting  virtual  patient  glucose  readings.  With  access  to  patient  data  from  Anne  Arundel 
Medical  Center,  this  combined  model  can  be  fit  to  real  patient  glucose  data  to  assess  parameter 
adaptation  success  and  verify  the  model’s  validity  for  use  as  the  foundation  of  a  personalized 
prescription  tool.  Pending  patient  data  confirmation,  this  combined  adaptive  model  is  prepared 
for  the  application  of  model  predictive  control  to  determine  appropriate  insulin  injection 
treatment  methods  for  diabetic  patients. 

The  topic  of  Diabetes  has  become  a  multigenerational  Trident  Research  focus  within  the  Systems 
Engineering  Department.  MIDN  Bryan  Weisberg  from  the  class  of  2013  investigated  the 
addition  of  glucagon  dynamics  to  the  Cobelli  model.  I  built  upon  his  knowledge  of  the  Cobelli 
model  to  incorporate  insulin  input  dynamics,  create  a  virtual  patient,  and  apply  parameter 
adaptation  to  fit  randomized  patient  data.  Next  year,  1/C  Alvin  Abes  will  apply  predictive  control 
to  this  adaptive  combined  model,  working  to  build  a  personalized  prescription  tool  which  can 
assist  physicians  in  treating  diabetic  patients.  These  continued  research  efforts  under  the 
instruction  of  Professor  Richard  O’Brien  will  further  the  development  of  a  control  algorithm 
which  will  improve  our  understanding  of  diabetes,  improve  patient  treatment  methods  with  the 
introduction  of  mathematical  prescription  tools,  and  contribute  to  the  development  of  an  artificial 
pancreas  in  the  future. 
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Enclosure  3:  Cobelli  Parameters 
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Units 

kl/kg 

1/min 

l/min 

1/min 

mg/kg 

l/kg 

1/min 

1/min 

1/min 

min*kg/pmol 

dimensionless 

dimensionless 

1/min 

1/min 

1/min 

1/min 

dimensionless 

mg/kg/m  in 

1/min 

pmol/l 

pmol/kg 

1/min 

Type  2 

Diabetic  Value 

(D 

tH 

0.042 

0.071 

0.0007 

cr> 

CD 

rsi 

o 

o 

0.379 

0.673 

0.269 

0.0526 

00 
t — i 
r—\ 
00 

o 

90 

0.0465 

CD 

o 

o 

o 

0.023 

0.0465 

6'0 

CD 

O 

o 

o 

O 

o 

00 

CD 

o 

0.00023 

60'0 

CT> 

o 

CO 

0.0007 

LO 

o 

o 

o 

0.0786 

9900'0 

Normal  Value 

00 

00 

T  1 

0.065 

0.079 

LO 

o 

o 

o 

CD 

CT> 

00 

00 

LO 

o 

o 

6T'0 

0.484 

0.194 

0.0304 

0.6471 

90 

00 

LO 

LO 

o 

o 

00 

o 

o 

o 

0.057 

00 

LO 

LO 

o 

o 

6'0 

0.00013 

CM 

00 

o 

0.00236 

t — 1 

O 

o 

2.7 

0.0021 

600'0 

00 

r—\ 

CD 

O 

O 

0.0079 

Physiological  Meaning 

distribution  volume  of  glucose 

glucose  rate  parameter  (b/n  plasma  and  tissue) 

glucose  rate  parameter  (b/n  plasma  and  tissue) 

glomerular  filtration  rate 

renal  threshold  of  glucose 

distribution  volume  of  insulin 

insulin  rate  parameter  (b/n  liver  and  plasma) 

insulin  rate  parameter  (b/n  liver  and  plasma) 

time  dependent  rate  parameter  (liver  degredation) 

insulin  rate  parameter  (periphery  degredation) 

Related  to  Hepatic  Extraction 

Related  to  Hepatic  Extraction 

Hepatic  extraction  of  insulin 

glucose) 

minimum  kemp 

rate  constant  of  intenstinal  absorption 

rate  of  grinding 

fraction  of  intestinal  absorption  which  appears  in  plasma 

percentage  of  dose  for  which  kemp  decreases  (flex  point 

percentage  of  dose  for  which  kemp  recovers  (flex  point  of 

EGP  at  zero  gluclose  and  insulin 

liver  glucose  effectiveness 

amplitude  of  insulin  action  on  liver 

amplitude  of  portal  insulin  action  on  liver 

delay  between  insulin  signal  and  action 
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mg/kg/m  in 

mg/kg/m  in 

pmol/l 

mg/kg 

l/min 

pmol/kg  per 

l/min 

pmol/kg/mirn  per 

l/min 

l/min 

rsj 

< 

3 

rsj 

< 
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< 
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o 
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O 
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transfer  rate  bewteen  portal  vein  and  liver 

trasform  rate  from  1  heameric  to  3  dimeric  molecules 

chemichal  equilbration  constant 

absorption  constant  parameter  (Lispro) 

shell  radius? 

insulin  degredation  rate  (Lipsro) 

Injection  amount 

rate  of  hexamer  disengagement  from  bound  state 

maximum  capacity  of  transformation 
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%Enclosure  4  -  Randomized  Virutal  Patient  Data 
%Creates  Random  patient  data  for  three  patients 

%  %show  original  cobelli  patient  data 
run  Diab_Parameter_Vector 

run  Parameter_initialization;  %picks  off  values  and  creates  variables 
sim  Cobelli_correct  %runs  Cobelli  Simulink  model 
figure  (1) 
plot (time, glue) 
clear  time  glue 


%change  each  parameter  by  a  random  percentage 


%set  acceptable  limits  for  the  different  parameters  (use  10%  of  difference 
%between  healthy  and  diabetic  in  cobelli  paper)  %corrected  to  be  centered 
%around  diabetic  value 


limit  =  [0.026174497 
0 . 054761905 
0 . 011267606 
0 . 028571429 
0 . 026022305 
0 . 025 

0 . 049868074 
0 . 02808321 
0 . 027881041 
0 . 042205323 
0 . 020288248 
0 

0.02 

0 . 005263158 
0 . 147826087 
0.02 
0 

0 . 116666667 
0 . 020588235 
0.002608696 
0.088888889 
0 . 012621359 
0.2 
0 . 08 

0 . 021374046 
0 . 01969697 
0 

0 . 046236559 
0 . 0265625 
0 . 051611935 
0 .060595238 
0 . 132323232 
0.284615385 
0 . 12 
0 
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0 

0 

] ; 


for  p  =  1:3 

change  =  zeros (37, 1) ; 
for  i  =  1:37 

change (i)  =  (2* (limit (i) )  . *rand (1, l)@limit (i) )  ;  %creates  a  vector  of  35  random 


end 


New_Param  =  zeros ( length (Parameters ),1); 


for  i  =  1:37 

New_Param(i)  =  ( 1+change ( i )) ^Parameters ( i ) ;  %modify  all  parameters,  using 

end 

Parameters  =  New_Param;  %reset  parameters  vector  to  randomized  parameters 

run  Parameter_init ializat ion;  %initialize  parameters  again  with  new  values 
sim  Cobelli_correct 


figure  (p+1) 
plot (time,  glue) ; 
xlabel ( ’ Time  (min) ' ) ; 

ylabel ( 1  Plasma  Glucose  Concentration  (mg/dl)'); 


if  p==l 

savefile  =  ' Pat ient l_3meals_correct . mat ' ; 
title ( 1  Plasma  Glucose:  Patient  1  Parameters'); 

end 


if  p==2 

savefile  =  ' Pat ient2_3meals_correct . mat ' ; 

t it le (' Plasma  Glucose:  Patient  2  Parameters'); 

end 


if  p==3 

savefile  =  ' Pat ient 3_3meals_correct . mat ' ; 

t it le (' Plasma  Glucose:  Patient  3  Parameters'); 

end 


save ( savefile, ' Parameters ' ,  'glue',  'time') 


%store  glucose  response 


end 


%plot  results 
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clear  all 
hold  on 

load  Pat ient l_3meals_Correct 
plot (time,  glue, ’r’) 
hold  on 

xlabel ( ’ Time  (min) ’ ) 

ylabel(’Gluc  (mg/dL) ’) 

clear  time  glue 

load  Pat ient2_3meals_Correct 

plot (time, glue, 'g') 

clear  time  glue 

load  Pat ient 3_3meals_Correct 

plot (time, glue,  ’b’) 

legend ( 1  Pat ient  lf,  ’Patient  2’,  ’Patient  3'); 

Error  using  Randomized_Data_Corrected  (line  7) 
Unable  to  load  block  diagram  ' Cobelli_correct  ’ 
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%  Enclosure  5  Part  1-  Parameter  sensitivity  testing  code 

%  changes  each  parameter  one  a  time,  runs  cobelli  model, 

%  and  calculates  root  mean  squared  error  in  the  glucose  output 
%  (a  similar  code  computes  and  saves  the  error  based  on  insulin  data) 

run  Parameter_vector ;  %  creates  normal  parameter  vector  and  defines  each  specific 

run  Parameter_init ializat ion;  %picks  off  values  and  creates  variables 
sim  Cobelli2  %runs  Cobelli  Simulink  model 
savefile  =  ' init_var . mat '  ; 

save ( savef ile,  ’glue',  'time')  %store  original  glucose  response 
init_gluc  =  glue; 


%change  each  parameter  one  at  a  time 

changes  =  [0.1,  0.2,  0.3,  0.4  ,0.5,  0.6,  0.7,  0.8,  0.9,  1];  %change  by  10%  increme 
results  =  zeros ( length ( ins ) ,  (( length ( changes ) *length (Parameters )) +1 )) ;  %create  sp 

results (:,1)  =  time; 

New_Param  =  zeros ( length (Parameters ), length ( changes )) ; 

Sens  =  zeros (( length (Parameters ) -2 ), length ( changes )) ; 
for  i  =  1:  ( length (Parameters ) -2 ) 
for  j  =  1 : length ( changes ) 
try 

New_Param ( i , j )  =  ( 1+changes ( j )) ^Parameters ( i ) ;  %modify  one  parameter  at  a 

%creates  a  matrix  where  each  row  represents  a  different  parameter 
%and  the  columns  show  the  changes  made 

run  Parameter_vector;  %reset  parameter  vector  from  previous  changes 
Parameters ( i )  =  New_Param ( i , j ) ;  %reset  parameter  vector  with  new  value 

run  Parameter_init ializat ion;  %initialize  parameters  again 
sim  Cobelli2  %run  Cobelli  model  with  new  parameters 

results  (: ,  (  (i-1) *10  + j  +  1) )  =  ins;  %create  matrix  of  glucose  results:  1st  c 

gluc_error  =  100*sqrt ( sum ( (gluc-init_gluc) . A2 ) / length (glue) ) / (max (glue) -mi 
Sens(i,j)  =  gluc_error/changes ( j ) ;  %similar  to  New_Param  -  shows  sensitivi 
i 

catch 

disp('Oh  No!  An  error  occured  -  execution  will  continue')  %  use  a  2  wh 
Sens  ( i , j )  =  2 ; 

end 

end 

end 


savefile  =  ' Results .mat ' ; 

save ( savef ile,  'results',  'Sens')  %save  glucose  results  in  matlab  file 

Error  using  run  (line  55) 

Parameter_vector  not  found. 

Error  in  Encl_4_Parameter_Sensitivity  (line  7) 

run  Parameter_vector ;  %  creates  normal  parameter  vector  and  defines  each  s 
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%Enclosure  5  Part  2:  Analyzing  Sensitivity  Results 

%combines  parameter  sensitivity  based  on  both  insulin  and  glucose  output  data 

%load  results  from  parameter  sensitivity  code 

load (  '  C : \Users\ml51 65  6\Document s\MATLAB\Sens_ins_full .mat  ' ) 

I  =  Sens; 

load (  '  C : \Users\ml51 65  6\Document s\MATLAB\Sens_gluc_full .mat  ' ) 

G  =  Sens; 

%  find  the  average  sensitivity  for  each  insulin  and  glucose 
maximum_I  =  zeros (35,  2); 
avg_I  =  zeros  (35, 2); 
for  i  =  1:35 

maximum_I ( i , 1 )  =  i; 

avg_I ( i , 1 )  =  i ; 

maximum_I ( i , 2 )  =  max ( I  ( i ,  :  )  )  ; 

avg_I ( i , 2 )  =  mean ( I  ( i ,  :  )  )  ; 

end 

maximum_G  =  zeros (35,  2); 
avg_G  =  zeros (35,2); 
for  i  =  1:35 

maximum_G ( i , 1 )  =  i; 

avg_G ( i , 1 )  =  i ; 

maximum_G ( i , 2 )  =  max(G(i,:)); 

avg_G(i,2)  =  mean (G ( i ,  : ) )  ; 

end 

%sort  results 

maximum_I ( : ,  3:4)  =  sort rows (maximum_I , -2 ) ;  %columns  3  and  4  will  show  sorted  valu 

avg_I ( : ,  3:4)  =  sortrows (avg_I, -2) ;  %columns  3  and  4  will  show  sorted  values 

maximum_G ( : ,  3:4)  =  sort rows (maximum_G, -2 ) ;  %columns  3  and  4  will  show  sorted  valu 

avg_G ( : ,  3:4)  =  sortrows (avg_G, -2) ;  %columns  3  and  4  will  show  sorted  values 

%combine  insulin  and  glucose  sensitivity 
Comb_max  =  zeros  (35,2); 

Comb_avg  =  zeros (35,2); 
for  i  =  1:35 

Comb_max ( i , 1 )  =  i; 

Comb_avg ( i ,  1 )  =  i; 

Comb_max ( i , 2 )  =  sqrt ( (maximum_G ( i , 2 ) ) A2+ (maximum_I ( i ,  2 ) ) A2 )  ; 

Comb_avg ( i , 2 )  =  sqrt ( (avg_G (i, 2 ) ) A2+ (avg_I (i,  2 )  )  A2 )  ; 

end 

%sort  results 

Comb_max ( : ,  3:4)  =  sort rows (Comb_max, -2 ) ;  %columns  3  and  4  will  show  sorted  values 
Comb_avg ( : ,  3:4)  =  sort rows (Comb_avg, -2 ) ;  %columns  3  and  4  will  show  sorted  values 

%save  in  excel  and  label  columns 
savefile  =  ' Sensitivity_Results_Gluc_Ins ' ; 

col_header  =  {’Parameter  Number’,  ’Max  Sensitivity ’,’ Parameter  Number’,  ’Max  Sensi 
row_header (1:35,1)  =  {'Vg',  ’kl',  ' k2 ' ,  ’kel’,  ’ke2',  ’VI’,  'ml',  'm2',  'm4',  'm5 


xlswrite ( savef ile,  Comb_max,  1,  'B2'); 

xlswrite ( savef ile,  col_header,  1,  fBlf); 
xlswrite ( savef ile ,  row_header,  1,  ’ A2 ’); 

col_header  =  {’Parameter  Number' ,  ’Avg  Sensitivity Parameter  Number' ,  'Avg  Sens 
xlswrite ( savef ile,  Comb_avg,  2,  'B2'); 

xlswrite ( savef ile ,  col_headerA  2,  'Bl'); 

xlswrite ( savef ile ,  row_header,  2,  ’A2'); 

Error  using  load 

Unable  to  read  file  ' C : \Users\ml51 65 6\Documents \MATLAB\Sens_ins_full .mat  ' 

Error  in  Encl_4_Sensitivity_Analysis  (line  5) 

load ( ’ C : \ Users \ml51 65  6\Documents \MATLAB\Sens_ins_full .mat  ' ) 
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%Enclosure  6  Part  1:  Parameter  Adaptation  Main  File 

%  This  is  the  top-level  script.  It  calls  the  optimization 
function  fmincon  and  makes  all  variable  declarations  necessary 
for  this  call.  Those  variable  declarations  include:  1)  The 
function  handle  for  the  "management"  function  (the  function  to 
be  called),  2)  The  upper  and  lower  bounds  for  each  parameter 
being  fit,  3)  Initial  condition  vector  for  the  parameters. 

These  elements  must  be  placed  in  the  appropriate  "spot"  in 
%  the  fmincon  function  call.  (Type  "help  fmincon"  in  Command 
Window  to  see  full  documentation  and  order  of  arguments.  Square 
braces  []  are  used  to  give  no  input  to  a  specific  function 
argument.  Order  matters,  and  fmincon  assigns  inputs  to  each 
argument  in  the  order  received,  so  use  []  to  skip  specific 
arguments  and  enter  later  ones.) 

%  Function  handle  for  function  to  be  called 
f  =  @ (p) Dilks_model_sim_10 (p) ; 

%  Lower  and  upper  bounds  of  parameters  (10  parameters  to  be 
fitted:  Vg, 

%  kl,  VI, m4,  m5,  m6,  kabs,  kgri,  kpl,  kp4)  (1,  2,  6,  9,  10,  11, 
15,  16,  22,  25) 

LB  =  [1.451  0.0397  0.039  0.2615  0.05038  0.79533  0.0196  0.04557 
3.051  0.07692]  ; 

UB  =  [1.529  0.0443  0.041  0.2765  0.05482  0.82827  0.0264  0.04743 
3.129  0.08028] ; 

%  Best  initial  guess 

pO  =  [1.49  0.042  0.04  0.269  0.0526  0.8118  0.023  0.0465  3.09 
0.0786] ; 

%  Calling  the  optimization  function 
[P  cost]  =  fmincon  (f , pO,  [ ] ,  [  ] ,  [  ] ,  [ ] , LB, UB) 


save  Parameter  Adaptation  Results  new  P  cost 
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%  Enclosure  6  Part  2:  Management  Function  (adapts  parameters) 

%  This  is  the  "management"  function.  It  handles  passing  the 
candidate 

%  parameter  vector  to  the  model  and  computing  the  cost  based  on 
the 

%  returned  model  simulation. 

o 

o 


function  F  =  Dilks_model_sim_10 (p) 

run  Diab_Parameter_Vector  %load  original  Diabetic  Cobelli 
parameter  values 

index  =  [1,  2,  6,  9,  10,  11,  15,  16,  22,  25];  %specific 
parameters  to  be  updated 
for  i  =  1 : length (index) 

Parameters (index (i) )  =  p(i);  %update  only  the  parameters  we're 

interested  in 

end 

run  Param_Adap_Setup  %assigns  new  parameter  values  and  creates 
new  Initial  conditions 

%uses  patient  data  to  create  time  vector  for  data  points 
Parameters (38)  =  lb; 

Parameters (39)  =  Sb; 

Parameters (40)  =  h; 

Parameters ( 4 1 )  =  meal; 

%  Printing  the  parameter  vector  chosen  by  fmincon  each 

interation 

P 

%  Solving  the  differential  equation  (running  the  model) 

[time,  estates]  = 

ode45 (@Cobelli, T_data, InitCond' , [ ] , Parameters) ; 
glue  =  estates (:, 9) /Vg; 


%  Calculate  the  cost  of  this  parameter  vector,  p,  by  comparing 
simulation 

%  to  data  at  the  same  time  points.  Print  this  value  each 
iteration . 

F  =  sum ( (glue-data) . A2 ) ; 


34 


%Enclosure  6  Part  3:  Setup  file  for  Parameter  Adaptation 

%Sets  up  parameter  variables,  initial  conditions,  and  data  time 

vector 

%Pick  off  Parameters  for  use  in  Simulink  Model 
Tstep  =  0.01; 

%  Glucose  Kinetics 
Vg  =  Parameters ( 1 )  ; 
kl  =  Parameters  ( 2 ) ; 
k2  =  Parameters  (3) ; 

%  Renal  Excretion 
kel  =  Parameters  ( 4 ) ; 
ke2  =  Parameters  ( 5 ) ; 

%  Insulin  Kinetics 

Vi  =  Parameters ( 6 ) ;  %note :  Vi  in  Weisberg 

ml  =  Parameters ( 7 ) ; 

m2  =  Parameters  (8) ; 

m4  =  Parameters  ( 9) ; 

m5  =  Parameters  ( 10 ) ; 

m6  =  Parameters  ( 1 1 ) ; 

HEb  =  Parameters ( 12 ) ;  %  where  is  this  used? 

%  Rate  of  Appearance  -  Gastro  Intestinal  Tract 

kmax  =  Parameters ( 1 3 ) ; 

kmin  =  Parameters ( 14 ) ; 

kabs  =  Parameters (15) ; 

kgri  =  Parameters (16) ; 

f  =  Parameters ( 17 )  ; 

a  =  Parameters  (18) ; 

b  =  Parameters (19) ; 

c  =  Parameters  (20) ; 

d  =  Parameters (21) ; 

%  Endogenous  Production  -  Liver 

kpl  =  Parameters (22 ) ; 

kp2  =  Parameters ( 2 3 ) ; 

kp3  =  Parameters ( 2 4 )  ; 

kp4  =  Parameters (25)  ; 

ki  =  Parameters ( 2 6 ) ; 

%  Utilization  -  muscle  and  adipose  tissue 
Fens  =  Parameters (27 ) ; 

VmO  =  Parameters (28) ; 

Vmx  =  Parameters (29) ; 

KmO  =  Parameters ( 30 ) ; 
p2U  =  Parameters (31) ; 

%  Secretion 
K  =  Parameters ( 32 )  ; 
alpha  =  Parameters (33) ; 

Beta  =  Parameters ( 34 ) ;  %note :  different  than  weisberg 
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gamma  =  Parameters (35) ;  %note :  different  than  weisberg 
%h  =  Parameters (36) ; 

%  Insensitive  Parameters 
%BW  =  Parameters ( 37 ) ; 

%D  =  Parameters ( 38 ) ; 

%  Prof  O'Brien  change  18  Nov 
BW  =  Parameters (36)  ; 

D  =  Parameters  (37); 


%Basal  Parameters  for  Cobelli  Model:  (initial  conditions) 

%  Insulin 
Sb  =  (m6  -  HEb)/m5; 
m3b  =  (HEb*ml) / (1-HEb) ; 
matrix  =  [ml+m3b,  -m2;  -ml,  m2+m4]; 
i_matrix  =  matrix^-l* [Sb;  0]; 
lib  =  i_matrix(l); 

Ipb  =  i_matrix(2); 
lb  =  Ipb/Vi ; 

Ipob  =  Sb/gamma; 

%  Glucose 
Ipob  =  Sb/gamma; 

syms  Gtb  Gpb  EGPb  %recognizes  unknown  variables  as  symbols 
%assuming  Gpb  >  Ke2 

51  =  ' kpl  ==  EGPb  +  kp2 *Gpb+kp3* Ib+kp4 * Ipob ' ;  %Equation  12 

52  =  'EGPb  ==  Fens  +  (VmO*Gtb)  /  (KmO+Gtb)  +  kel*  (Gpb  -  ke2)'; 

%Equation  14,  15,  27 

53  =  ' (VmO*Gtb) / (KmO+Gtb)  =  kl*Gpb  -  k2*Gtb';  %Equation  1 
S_Gluc  =  solve (SI,  S2,  S3,  Gtb,  Gpb,  EGPb);  %returns  a 

structure  with  solved  equations 

tmp  =  double (subs (S_Gluc . Gtb) ) ;  %substitute  in  known  parameters 
Gtb  =  tmp (2);  %  pick  out  positive  solution 
tmp  =  double (subs (S_Gluc . Gpb) ) ; 

Gpb  =  tmp ( 2 ) ; 

tmp  =  double ( subs (S_Gluc . EGPb) ) ; 

EGPb  =  tmp ( 2 ) ; 

if  (Gpb<=ke2)  %checks  if  Gpb  <  ke2  -  if  so,  need  to  change 
equations 

clear  Gtb  Gpb  EGPb; 
syms  Gtb  Gpb  EGPb; 

51  =  'kpl  ==  EGPb  +  kp2*Gpb+kp3*Ib+kp4*Ipob ' ;  %Equation  12 

52  =  'EGPb  ==  Fens  +  (VmO*Gtb) / (KmO+Gtb) ' ;  %Equation  14, 

15,  27 

53  =  ' (VmO*Gtb) / (KmO+Gtb)  =  kl*Gpb  -  k2*Gtb';  %Equation  1 
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S_Gluc  =  solve (SI,  S2,  S3,  Gtb,  Gpb,  EGPb) ;  %returns  a 
structure  with  solved  equations 

tmp  =  double (subs (S_Gluc . Gtb) ) ;  %substitute  in  known 
parameters 

Gtb  =  tmp (2);  %  pick  out  positive  solution 
tmp  =  double (subs (S_Gluc . Gpb) ) ; 

Gpb  =  tmp ( 2 ) ; 

tmp  =  double ( subs (S_Gluc . EGPb) ) ; 

EGPb  =  tmp ( 2 ) ; 

end 


%  Prof  O'Brien  edit  18  Nov 
Gb  =  Gpb/Vg; 
h  =  Gb; 


%  meal  at  t=l 
meal  =  1; 


zeros  (12,1);  % 


InitCond  = 
InitCond ( 1 ) 
InitCond ( 2 ) 
InitCond ( 3 ) 
InitCond  ( 4 ) 
InitCond  ( 5 ) 
InitCond (6) 
InitCond ( 7 ) 
InitCond ( 8 ) 
InitCond  (9) 
InitCond (10) 
InitCond ( 11 ) 
InitCond  ( 12 ) 


=  D;  %Qstol 
=  0;  %Qsto2 
=  0;  %Qgut 
=  lb;  %I1 
=  lb;  %Id 
=  Ipb;  %Ip 
=  Ipob;  %Ipo 
=  0 ;  %X 
=  Gpb;  %Gp 
=  Gtb;  %Gt 
=  0;  %Y 
=  I lb;  %I1 


initial  values 


-  find  this ! 


for  all  states 


load  Patientl_correct_data  %  load  random  patient  data 
%record  what  the  original  patient  parameters  were  for 
verification 

%original  parameters:  1.4636,  0.0408,  0.0394,  3.0783,  0.0798 
%  pick  off  specific  time  points  from  the  data  -  glucose  readings 
every  20 

%  min  (15  readings  total) 

T_data  =  zeros  (16,1); 
data  =  zeros (16,1); 

T_data(l)  =  pat_time(l); 
data(l)  =  pat_gluc(l); 
for  i  =  1:15 

T_data(i+1)  =  pat_time (i*20*1000) ; 
data(i+l)  =  pat_gluc (i*20*1000) ; 


end 
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%Enclosure  6  Part  4:  Cobelli  Model  Equations 
%runs  the  cobelli  model  using  an  ODE  solver 


%inputs : 

%t=current  time 

%state=current  state  of  tracked  variables 
%Parameters=parameter  vector 

%inputs  =  2  row  vector  of  insulin  and  glucagon  injection  amounts 


function  stateReturn  =  Cobelli (t,  state.  Parameters) 


%%  plucks  off  parameters,  assigns  variable  names,  creates 
initial  conditions 

%Pick  off  Parameters  for  use  in  Cobelli  model 

%  Glucose  Kinetics 
Vg  =  Parameters (1) ; 
kl  =  Parameters  (2 ) ; 
k2  =  Parameters  (3) ; 

%  Renal  Excretion 
kel  =  Parameters ( 4 ) ; 
ke2  =  Parameters (5) ; 

%  Insulin  Kinetics 
Vi  =  Parameters ( 6) ; 
ml  =  Parameters ( 7 ) ; 
m2  =  Parameters  (8) ; 
m4  =  Parameters  ( 9) ; 
m5  =  Parameters  (10) ; 
m6  =  Parameters  (11) ; 

HEb  =  Parameters  ( 12 ) ; 

%  Rate  of  Appearance  -  Gastro  Intestinal  Tract 
kmax  =  Parameters (13) ; 
kmin  =  Parameters ( 14 ) ; 
kabs  =  Parameters (15) ; 
kgri  =  Parameters (16) ; 
f  =  Parameters ( 1 7 )  ; 
a  =  Parameters ( 18 )  ; 
b  =  Parameters (19)  ; 
c  =  Parameters (20)  ; 
d  =  Parameters (21)  ; 

%  Endogenous  Production  - 
kpl  =  Parameters (22 )  ; 
kp2  =  Parameters (23)  ; 


Liver 


kp3  =  Parameters ( 2 4 ) ; 
kp4  =  Parameters (25)  ; 
ki  =  Parameters (2  6)  ; 

%  Utilization  -  muscle  and  adipose  tissue 
Fens  =  Parameters (27 )  ; 

VmO  =  Parameters (28) ; 

Vmx  =  Parameters ( 2 9 ) ; 

KmO  =  Parameters ( 30 ) ; 
p2u  =  Parameters  (31) ; 

%  Secretion 
K  =  Parameters ( 32 ) ; 
alpha  =  Parameters (33) ; 

Beta  =  Parameters ( 34 )  ; 
gamma  =  Parameters (35) ; 

%  Prof  O'Brien  change  18  Nov 
BW  =  Parameters (36) ; 

D  =  Parameters ( 37 ) ; 
lb  =  Parameters (38) ; 

Sb  =  Parameters (39) ; 
h  =  Parameters (40)  ; 
meal  =  Parameters (41)  ; 

%%  State  Vector 
%Rate  of  Appearance 

Qstol  =  state (1);  %where  meal  input  is  implemented  as  initial 

condition 

Qsto2  =  state (2); 

Qgut  =  state (3); 

II  =  state  ( 4 )  ; 

Id  =  state  (5);  %delayed  insulin  signal 
Ip  =  state ( 6 ) ; %plasma  insulin  concentration 
Ipo  =  state ( 7 ) ; 

X  =  state (8);  %interstitial  insulin  (muscle  and  adipose) 

Gp  =  state ( 9)  ; 

Gt  =  state ( 10 ) ; 

Y  =  state (11) ; 

II  =  state  ( 12 )  ; 


%%  Model  Equations  %does  the  order  of  these  equations  matter 
%Rate  of  Appearance 
Qsto=Qstol+Qsto2 ; 

KemptQsto=kmin+ (kmax-kmin) /2* (tanh (a* (Qsto-b*D) ) -tanh (c* (Qsto 
d*D) ) +2 )  ; 

dQstol=-kgri*Qstol ; 

dQsto2=-KemptQsto*Qsto2+kgri*Qstol ; 
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dQgut=-kabs*Qgut+KemptQsto*Qsto2 ; 
Ra=f *kabs*Qgut/BW; 

%Liver 

I=Ip/Vi; 

dll=-ki* (Il-I) ; 

dld=-ki* (Id-11) ; 

EGP=kpl-kp2*Gp-kp3*Id-kp4*Ipo; 

%Muscle  and  Adipose  Tissue 
dX=-p2u*X+p2u* (I-Ib) ; 
VmX=VmO+Vmx*X; 

Uid=VmX*Gt/ (KmO+Gt ) ; 

Uii=Fcns ; 


%Glucose  Kinetics 
if  (Gp>ke2) 

E=kel* (Gp-ke2)  ; 

else 

E=0 ; 

end 

dGp=EGP+Ra-Uii-E-kl*Gp+k2*Gt; 

dGt=-Uid+kl*Gp-k2*Gt; 

G=Gp/Vg / 
dG=dGp/Vg; 

%Pancreas 

%static  secretion 
if (Beta* (G-h) >=-Sb) 

dY=-alpha* (Y-Beta* (G-h) ) ; 

else 

dY=-alpha*Y-a*Sb; 

end 


%dynamic  secretion 
if (dG>0 ) 

Spo=Y+K*dG+Sb; 

else 

Spo=Y+Sb; 

end 

dIpo=-gamma*Ipo+Spo; 

S=gamma*Ipo; 

%Insulin  Dynapics 
HE=-m5*S+m6; 
m3=HE*ml/ (1-HE) ; 


dll=- (ml+m3) *Il+m2*Ip+S; 
dlp=- (m2+m4) *Ip+ml*Il; 


%%  return 

stateReturn= [dQstol ; dQsto2 ; dQgut; dll ; did; dip; dlpo; dX; dGp; dGt; dY 

dll]; 

end 

%this  function  returns  state  derivatives  and  uses  them  and 
the  current 

%state  to  determine  the  next  state 


